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Nonlinear quantum optical systems are of paramount relevance for modern quantum technolo¬ 
gies, as well as for the study of dissipative phase transitions. Their nonlinear nature makes their 
theoretical study very challenging and hence they have always served as great motivation to develop 
new techniques for the analysis of open quantum systems. In this article we apply the recently 
developed self-consistent projection operator theory to the degenerate optical parametric oscilla¬ 
tor to exemplify its general applicability to quantum optical systems. We show that this theory 
provides an efficient method to calculate the full quantum state of each mode with high degree of 
accuracy, even at the critical point. It is equally successful in describing both the stationary limit 
and the dynamics, including regions of the parameter space where the numerical integration of the 
full problem is significantly less efficient. We further develop a Gaussian approach consistent with 
our theory, which yields sensibly better results than the previous Gaussian methods developed for 
this system, most notably standard linearization techniques. 

PACS numbers: 42.50.-p, 03.65.Yz, 42.65.Sf, 42.65.Yj 


I. INTRODUCTION 

Nonlinear optical systems play an important role in 
the field of optics both in classical BE] and in quantum 
(3H3 regimes. Quantum mechanical effects, in particu¬ 
lar, which are not explainable by classical optics, have 
triggered substantial research, especially in connection 
to modern applications such as high-precission measure¬ 
ments jBHEE] and quantum information communication 
and processing PJEE3]. Importantly, the nonlinear na¬ 
ture of these systems leads to non-Gaussian states, which 
typically precludes an analytic treatment and therefore 
requires elaborate theoretical approaches mm- 

In a system where the dynamical degrees of freedom 
evolve on different time scales, approximate descriptions 
of reduced complexity may be found. For example, adi¬ 
abatic elimination techniques can be exploited to derive 
effective equations of motion mm- In this work, we 
apply the recently introduced self-consistent projection 
operator theory HU to the degenerate optical paramet¬ 
ric oscillator , and exemplify how it generalizes adiabatic 
elimination approaches. This theory takes dynamical 
back-action between the degrees of freedom into account 
and therefore does not require any time-scale separation. 
We expect our method to be directly applicable to other 
nonlinear quantum optical models such as those for non¬ 
degenerate or multi-mode parametric oscillation fT9H22j . 
lasing m m na mj, optomechanical parametric oscilla¬ 
tion [25j [26] , or the dissipative Dicke model [271129] . 
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Degenerate optical parametric oscillators (DOPOs) 
have been extensively studied in the past mm and are 
one of the paradigm examples of a system subject to a 
driven and dissipative phase transition. It is formulated 
as a bosonic problem with two modes, signal and pump , 
subject to dissipation and interacting nonlinearly. In the 
adiabatic limit of a fast decaying pump mode, an effec¬ 
tive master equation can be derived by means of standard 
projection operator approaches [15] and due to its re¬ 
duced complexity, the steady state can be found by solv¬ 
ing the corresponding Fokker-Planck equations for the 
positive P distribution nano. Yet away from the adia¬ 
batic limit one has to resort to numerical simulations or 
perturbative treatments [321137] . Non-equilibrium many- 
body techniques such as the Keldysh formalism have also 
been employed to study steady-state properties [38ti40] . 
While the application of all these techniques has allowed 
to deepen our understanding of DOPOs and phase transi¬ 
tions in driven dissipative quantum systems enormously, 
it is important to note that they are naturally built to 
determine the evolution of observables, making the deter¬ 
mination of the quantum state of the optical fields very 
challenging, if not impossible. 

Our approach, in contrast, derives a set of coupled 
equations for the reduced states of the two optical modes 
of the DOPO. By numerically solving these equations, we 
find the reduced density matrices of both the pump and 
the signal modes. We test the accuracy of our method by 
comparing its results with those of the full DOPO prob¬ 
lem in regions of the parameter space where this is nu¬ 
merically tractable. Our findings show that our method 
is remarkably close to the exact results, both for steady 
states and dynamics, while being less numerically de¬ 
manding than the full simulation of the DOPO problem. 
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FIG. 1. Sketch of the self-consistent projection operator the¬ 
ory for the DOPO, which consists of an optical cavity contain¬ 
ing a crystal with second-order optical nonlinearity, pumped 
by a laser at frequency w p — 2 uj s (pump mode), and capable 
of producing a field at the subharmonic frequency uj s (signal 
mode) via down-conversion in the crystal. For the sake of 
illustration, we consider in the figure a doubly-resonant semi- 
monolithic configuration in which each face of the crystal acts 
as a mirror for one of the modes, but is transparent for the 
other, allowing to create independent cavities for the pump 
and signal modes via two additional partially transmitting 
mirrors 46 . In our approach, the full problem described by 
the state p(t) and the Liouvillian C is mapped onto two cou¬ 
pled equations for the signal and pump modes. In one of the 
equations, the signal mode considered as the system is de¬ 
scribed by an effective master equation for its reduced state 
p s (t) = £ s (pp)p s (t) with an effective Liouvillian depending on 
the state of the pump, which plays here the role of an envi¬ 
ronment. The other equation considers the reversed scenario 
with the pump taking the role of the system while the sig¬ 
nal is interpreted as the environment leading to the effective 
equation p p {t) = £p(p s )p P (t). In this way the two equations 
form a closed set. 


qualitative picture of the physics in many, albeit not all, 
systems, it leads to unphysical predictions close to the 
critical points of the classical theory, e.g., to infinite pho¬ 
ton numbers in the case of the DOPO [42]. These unphys¬ 
ical predictions can be regularized by applying a more 
elaborate Gaussian state approximation where the sys¬ 
tem is not forced to stay in its classical state, but chooses 
instead an average configuration more consistent with the 
quantum fluctuations that perturb it [43]. Motivated by 
such an idea, we apply a Gaussian approximation within 
the self-consistent projection operator theory, and show 
that it gives more accurate quantitative results than any 
of the usual Gaussian techniques, as it does not assume 
a Gaussian state for the entire system, but only for the 
reduced state of one of the modes. 


The remainder of the paper is organized as follows. 
In Sec. un we introduce the DOPO model. We also dis¬ 
cuss its symmetries and briefly elaborate on the stan¬ 
dard linearization approach in Sec. IIA Sec. Ill reviews 


the main concepts of the self-consistent projection oper¬ 
ator theory and introduces the self-consistent Mori pro¬ 
jector (c-MoP) equations, which lie at the center of our 
study. Our theory provides a systematic extension of 
mean-field approaches as demonstrated in Sec. |III A] and 
reproduces known results in the adiabatic and the dia- 
batic limits introduced in Sec. Ill C An efficient proce¬ 
dure designed to deal with the non-Markovian structure 
of the c-MoP equations is provided in Sec. |IIIB[ which 
we use in Sec. [lV]to test the accuracy of our method for 
steady-state quantities and to present quantum states of 
the signal mode. A Gaussian state approximation on the 
c-MoP equations is performed in Sec. [V] which is shown 
to lead to highly accurate quantitative results as com¬ 
pared to previous linearization techniques. As a further 
test, we check in Sec. VI that our method provides the 
same level of accuracy for the dynamics, as it does for 
steady states. Finally, we conclude our work and present 
an outlook in Sec. Em 


It thus gives access to the reduced states of the cavity 
modes in regions of the parameters that are inaccessible 
to the latter. 

The possibly largest reduction of complexity in nonlin¬ 
ear quantum optical systems, however, comes from the 
application of Gaussian approximations on the state of 
the system. Within a Gaussian theory one can basically 
cover the whole parameter space efficiently to determine 
both steady-state and dynamical quantities such as two- 
time correlation functions. The simplest and most widely 
used Gaussian approach is known as the linearization 
technique [301111], which consists in assuming that the 
system configuration is, on average, in its classical state, 
but is constantly driven out of it by some “small” quan¬ 
tum fluctuations. While this technique provides a good 


II. THE DEGENERATE OPTICAL 
PARAMETRIC OSCILLATOR 

A DOPO consists of a driven optical cavity contain¬ 
ing a crystal with second order optical nonlinearity, see 
Fig. 0 Two relevant resonances at frequencies uj s ( sig¬ 
nal mode) and uj p = 2u) s (pump mode) exist in the cav¬ 
ity, which are nonlinearly coupled via parametric down- 
conversion inside the crystal, capable of transforming 
a pump photon into a pair of signal photons, and vice 
versa. We assume that the external driving laser is res¬ 
onant with the pump mode. Including damping through 
the partially transmitting mirrors at rates and for 
the pump and signal modes, respectively, the equation 
governing the evolution of the state p of the system in a 
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picture rotating at the laser frequency is given by pacng 

p(t) = [ e p( a l - a p) + |(<vT - ala 2 s ),p(t) 

+ 7 ? [ 2 a,jp(t)a] - a]a jP (t) - p{t)a ] aj \, ^ 

3=3,P 

where y /2 is the down-conversion rate and e p is propor¬ 
tional to square root of the injected laser’s power. We 
have defined bosonic operators a p and a s for the pump 
and signal modes, respectively, which satisfy canonical 
commutation relations [aj, a l} = s ji and [a,j. a{\ = 0 . 
Note that the nonlinear interaction is third order in the 
field operators, precluding a general analytic solution of 
Eq. § to which we refer as the Liouville-von Neumann 
equation or simply the full master equation of the DOPO. 


A. Linearization approach and symmetry breaking 

The right hand side of Eq. 0 can also be written in 
a shorthand notation by introducing a superoperator C 
(Liouvillian ), such that p(t) = Cp(t). For the major part 
of this work, we will be interested in the steady state 
p ss = lirn^oo p(£), which fullfills the equation Cp ss = 0. 
Due to the dissipation acting on both modes and because 
an arbitrarily large but finite truncation will always pro¬ 
vide an arbitrarily good approximation, we expect the 
steady state to be unique mm • 

We further note the invariance of the Liouvillian un¬ 
der a unitary transformation U 2 of Ising-type Z 2 which 
transforms a s as U 20 l s U\ = — a s . Since the steady state 
is unique, this implies that it has to be invariant under 
the Z 2 transformation too, i.e. U 2 PSSU 2 = Pss • This 
in turn leads to vanishing steady state expectation val¬ 
ues which include odd powers of the signal field operator 
a s . In particular (a s ) = 0 = (a p a\), as for example 
(a s ) = Tr{a s / 9 ss } = Tr{£7 2 as = ~( a s)- 

However, the most common technique used to analyze 
Eq. Q, known as the linearization approach , breaks this 
Z 2 symmetry [301I4T] . which has to be restored “by hand” 
at the end of the calculation, following the procedure 
that we explain at the end of Sec. |IV| Even though this 
method is more naturally introduced in the Heisenberg 
picture using the language of quantum Langevin equa¬ 
tions, it also admits a Schrodinger picture interpretation 
in terms of two successive approximations in the master 
equation. It starts by writing the bosonic operators as 
a j = a j + 5a j, with a.j = (aj) and hence (Saj) = 0 . 
In the first approximation, the fluctuation operators 5a j 
are neglected altogether; the evolution equations for (aj) 
(Bloch equations), then provide a set of nonlinear differ¬ 
ential equations for the amplitudes oq, which in the case 
of the DOPO read 

. _ X 2 

OL p — € p 7 P a p 2 ( 2 ) 

d s = + XOLpOLg. 


These correspond to the classical equations of the sys¬ 
tem, as they could have been obtained directly from 
Eq. 0 by assuming a coherent state for p ss , or simply 
from Maxwell’s equations. Depending on the injection 
parameter a = X e p/Xsjp one finds two types of steady- 
state solutions of Eq. ([2]). One of them has a s = 0 and 
ol v = e p /x P , and hence it does not break the symmetry; 
it is known as the below-threshold solution , and is only 
stable for a < 1. The o ther solution is bistable and has 
X a s = ±a/ 2 (xe p - 7 s) and %a p = j 3 , hence breaking 
the Z 2 symmetry; it is known as the above-threshold so¬ 
lution , and exists only for a > 1. The threshold point 
o = l marks a critical point where the classical theory 
predicts a phase transition from a signal-off phase with 
a s = 0 to a signal-on phase with a s 7 ^ 0. In the signal- 
off phase all injected power e p goes into the pump mode, 
while after crossing the critical point all the extra injec¬ 
tion is transferred to the signal mode, see the gray thin 
solid line in Fig. [2] 

Once the classical solutions have been identified, the 
second approximation consists in coming back to the orig¬ 
inal master equation with the bosonic operators written 
as aj = Oij + 5aj, and neglecting any term which goes 
beyond quadratic order in the fluctuation operators 5aj. 
This leads to a so-called linearized master equation which 
can be easily solved. 

One has to keep in mind that this linearized theory 
can only be trustworthy when the classical solution is a 
strong attractor, because only then the quantum fluctu¬ 
ations driving the system out of equilibrium are strongly 
damped, and quantum noise can be treated as a small 
perturbation. This means that, in particular, any pre¬ 
dictions obtained through this method cannot be trusted 
in the vicinities of critical points of the classical theory: 
points of the parameter space where one solution becomes 
unstable, making way for a new solution to kick in, hence 
creating non-analytic behaviour in some observable, that 
is, a classical phase transition. Indeed, this is exactly the 
case for the DOPO, in which this linearized description 
breaks down at threshold, offering unphysical predictions 
such as infinite photon numbers in the signal field (as il¬ 
lustrated by the gray thin line in Fig. [4]). 


III. SELF-CONSISTENT MORI PROJECTOR 
APPROACH 

To explain the approach employed in our calculations, 
we will first recapitulate some basic ideas of the self- 
consistent projection operator theory m- The first step 
is to divide the entire system into subsystems. In the 
DOPO this naturally amounts to consider the pump 
mode described by its reduced state p p (t ) = Tr s {p(t)} 
and the signal mode described by p s (t) = Tr P {p(t)}. In 
the spirit of open system theory jH EH we will first treat 
the pump mode as an “environment” for the signal mode, 
which then takes the role of the open “system”. Techni¬ 
cally this is done by introducing the time-dependent, self- 
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consistent Mori projector 'Pf(-) = p p (t ) ( 8 ) Tr p {-} whose 
action on the full state p(t) gives the factorized state 
Vfp(t) = p p (t) ( 8 ) p s (t). The term “self-consistent” is 
chosen because the state of the pump in Vf is not a 
time-independent reference state but is rather obtained 
consistently from the time-evolving state p(t) of the full 
dynamics. Using this projector, we derive a generalized 
Nakajima-Zwanzig equation which is an exact equation 
for the reduced state of the signal mode [18] . The ef¬ 
fective Liouvillian describing such a Nakajima-Zwanzig 
equation will depend on the state of the pump p p (t). In 
order to obtain a closed set of equations we need to re¬ 
verse the scenario and treat the pump mode as the “sys¬ 
tem” and the signal mode as the environment, see Fig. [l] 
for an illustration. 

Again analogous to open system theory, we split the 
full Liouvillian C from Eq. 0 into three parts. After 
performing a displacement a p —>> a p + ol p , where a p will 
be chosen later, see Sec. Ill A[ we write C = C p + C s + £/, 
with 


£p(’) — 7 P & P ) a P ( e P , • ] +7 P D ap (’) 

£ s (-) = ^M 2 -«;«s,-]+7A(-) 

£/(•) = | K al 2 -ala 2 s , • ] , 

(3) 


where we have defined the standard Lindblad superoper¬ 
ator Pfc(-) = 2-b^b(-) — (-)tfb, with b being an arbi¬ 
trary operator. The displacement a p a p -\-a p moves the 
large coherent background of the pump field into the free 
evolution of the signal C S: keeping only the pump mode’s 
fluctuations within the nonlinear signal-pump interaction 
£/. Such a step is important as our theory expands in 
powers of the interaction Liouvillian £/ in order to solve 
the Nakajima-Zwanzig equation. As in reference m we 
will expand to second order in the system-environment 
interaction. This approximation is known as the Born 
approximation [24] . The effective equations of the signal 
and the pump mode then read, 


Ps(t) = c s ps{t) +1 [d 2 (%>)(*) - “ 2 (%>)*(*)> a* (*)] ( 4 ) 

pt 


+ (|) | a 2 s ,J dt'e Cs{t tf) lC s (t,t f )p s (t') 


H.c. 


Pp(t) = £ P Pp(t) + | [a P (a 2 J(t) - a ] p {a 2 s )(t ), p p {t)\ (5) 


(f) { Up, f dt ' eCp{t 


+ £ 


+ H.c. 


where we have defined the Kernel superoperators 


= 8a p (t')(-)df(t,t') - (-)Sa p (tj d+ (t,t') (7) 
- 8al(t ')(•) d~(t, tj + (-)Sal(t') d~(t, t'), 


and for any operator Aj acting on the signal (j = s) or 
pump (j = s) subspace, we have defined the correspond¬ 
ing fluctuation operator SAj(t ) = Aj — Tr j{Aj Pj(t)}. 

The state of the pump mode p p (t) enters the signal 
mode’s dynamics, eq. Q, via ( a p )(t ) = Tr P {a p p p (t)} and 
the correlation functions 


d+(t,tj = Tr p { ale c ^-^8al(tjp p (t')}, (8) 

d p (t,t') = Tr p {ale Cp( ' t ~ t ' ) p p (tj8al(tj}, 

= Tr p {a i p e c ^ t ~ t '' , 5a p (tjp p (tj}, 
dp{t,tj = Tr p {ale c ^ t ~ t ' ) p p (tj6a p (tj}, 

In turn, the state of the signal mode p s {t) enters the 
pump mode’s dynamics, eq. 0 , via the expectation value 
(< a 2 s )(t ) = Tr s {alp s (t)} and the correlation functions 


= Ttsial^-OSal^psit')}, (9) 
df(t,tj = T8c s {al 2 e Cs{ - t ~ t ') Ps (tjba\ 2 (tj}, 
d~(t,tj = Tr s {a\ 2 e Cs(t - t ')5a 2 s (tjp s (tj}, 
d~(t,tj = Tr s {al 2 e c ^ t ~ t ' ) p s (tj5a 2 s (tj}. 


Equations 0 and 0 should be understood as two 
coupled equations which represent effective equations for 
the reduced states of the signal and the pump mode. We 
refer to these two equations as the c-MoP (consistent 
Mori Projector) equations of the DOPO. They can be 
thought of as non-Markovian and nonlinear master equa¬ 
tions which do not rely on any time-scale separation be¬ 
tween the modes. We will elaborate in detail on the limits 


where time-scale separation is present in Sec. Ill C 


The only assumptions made so far are the Born ap¬ 
proximation and the assumption of an initially factor¬ 
ized state p(0) = p p ( 0) (8) p s ( 0). The latter seems very 
reasonable by considering the vacuum as the state of the 
two modes before the driving laser is switched on. We 
also emphasize, our approach does not ignore system- 
environment or rather signal-pump correlations. In fact, 
it has been shown [18] that the Born term , the term sec¬ 
ond order in £/ which is here proportional to (y/2) 2 , 
clearly takes signal-pump correlations into account. We 
will show the crucial importance of the Born term in sev¬ 
eral examples below. Of course, c-MoP theory or any the¬ 
ory based on the concept of projection operators does not 
give access to explicit expressions for system-environment 
correlation functions. An example in this context could 
be the cross-correlation function (aja s ) — (aj,)(a s ). 

The most striking advantage of projection operator 
theories and in particular of the c-MoP theory is the 
reduction of the complexity of the problem. In the ex¬ 
ample of the DOPO the complexity of the Liouville-von 
Neumann eq. M scales as dim H s x dimH p , where H s / P 
denotes the Hubert space of the signal/pump modes, 
while the complexity of the c-MoP equations scale as 
dim% s + dim ji p . The self-consistent Mori-projector the¬ 
ory thus offers a very significant reduction of complexity. 
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A. Mean-field Approximation 


B. Born terms 


A merely approximate but very simple way of solving 
the c-MoP equations is to consider all terms up to first 
order in the interaction Cj only. Hence we drop all terms 
proportional to y 2 from eqs. Q and Within this 
approximation it does not make a difference whether the 
pump field is displaced or not. For simplicity we put 
the displacement a p from eq. ^ to zero and obtain two 
coupled equations 


pp(t) 

Ps(t) 


( e p 2 ) a p H.c. , p p {t) +7 p'DdpPpif)’* 


X 


[( a p )al 2 - H.c. ,p a (t)] +7 s^a s p s (t), 


( 10 ) 


known as mean-field equations [48]. These equations 
are quadratic in the field operators and therefore it is 
straightforward to solve them either numerically for the 
dynamics or analytically for the fixed points [43j[48]. The 
stationary state of the signal mode will be a Gaussian 
state H2Q3HH] centered around a vanishing field ampli¬ 
tude ( a s ) =0 as the mean-field equations do not break 
the Ising-type Z 2 symmetry. The steady state of the 
pump mode will be a coherent state with an amplitude 
given by (a,)?/= (e,-$(a*)?/)h P . 

Just like the c-MoP equations Q and <§ , the mean- 
field equations are coupled nonlinear equations which 
have to be solved self-consistently. Within mean-field 
theory fluctuations of the pump mode are disregarded. 
Fluctuations of the signal mode, however, are (at least 
to some extend) taken into account [43, 48]. This leads 
to the regularization of the divergences appearing in the 
classical theory or rather the standard linearization ap¬ 
proach. For our purposes it is important to note that 
the pump field amplitude always stays below the clas¬ 
sical above-threshold solution, i.e. {a p )^ F < 7 s /x- I n 
the remainder of the paper we will use it as the displace¬ 
ment in eq. O), i.e. a p = {a p )^ F . This will guarantee 
a well-behaved Liouvillian for the free system C s as we 


will explain in more detail in Sec III B 


The mean-field equations can also be found by putting 
the factorized state Ansatz p(t) = p p {t) (8) p s (t) into the 
Liouville-von Neumann equation, here given by eq. 
before tracing out each of the modes separately. This 
well-known procedure, indeed, neglects all signal-pump 
correlations. Within the self-consistent projection op¬ 
erator theory, mean-field can be understood as an ap¬ 
proximation to linear order in the interaction Cj for the 
dynamics of reduced density matrices. Our theory there¬ 
fore provides a systematic generalization of mean-field 
approaches. It is due to the Born terms, which are sec¬ 
ond order in £/, that signal-pump correlations are taken 
into account. Therefore, we expect a different quality of 
approximation by going from first order to second order 
in the interaction. 


In order to solve the full c-MoP equations including the 
Born terms we will need to overcome two main difficul¬ 
ties. While the c-MoP equation of the pump mode is 
quadratic in the field operators, granting us with a closed 
set of equations including only first and second moments 
of the pump field, the c-MoP equation of the signal 
mode is quartic in the field operators. We will there¬ 
fore either solve the equation of the signal fully numeri¬ 
cally, see Sec. El or apply a Gaussian state approxima¬ 
tion as presented in section Sec. [VJ In any of these two 
approaches, we need to overcome the second difficulty 
which arrises due to the non-Markovian structure of our 
theory. In the remainder of this section we will thus show 
how to rewrite an integro-differential equation of first or¬ 
der into a set of coupled ordinary differential equations. 
For the present problem this step is crucial, as solving 
the integro-differential equations is significantly more de¬ 
manding for both numerical and analytical approaches. 

We start by evaluating the correlation functions of the 
pump. By taking derivatives of the pump correlators 
dp (£, t') and d p (£, t') with respect to £, see eq. |8|, consid¬ 
ering initial conditions at t = t' (note that we understand 
from the c-MoP equations that t' < £), and exploiting the 
fact that the operator Sa^ p (t')p p (t') is traceless, we find 

dp(t,t') = dp (t, t') = [(<4 2 )(i') - <a P )* 2 (0] e-7p(t-t,) > 

d-(tx) = [1 + Ha P m - ikxoiv ^-7 

d-(t,f) = Ha p )(t ') - IKHOIV 7 ^-* 0 - 

( 11 ) 

Hence, all correlation functions of the pump can be writ¬ 
ten in a form where the t dependence only enters in a 
simple exponential factor. 

A bit more effort is needed in order to simplify the 
correlation functions of the signal, but the main steps 
are mainly identical. All the functions in eq. are of 
the form /(t, t') = Tr s {a\ 2 e Cs( d~ t ^ A{t')} with a traceless 
operator A(t r ) depending solely on t'. Again, we take the 
derivative of f(t, t') with respect to t and find an equation 
of motion of the form dtVt'(t) = Mv t '(t) with a column 
vector 


Vf{t) = col ((aja s ), (a 2 ), (aj 2 )^) , (12) 

where the expectation values with the tilde are defined in 
the usual way as the trace over the signal mode but with 
a density matrix given by pt'(t) = ^A(P). The 

matrix M reads 


f -2q s xa P X&l \ 

M = 2 X a* p —2 7s 0 

\ 2 xa p 0 - 27 S J 

It is straight forward to diagonalize M. We write 
M = I/AI/ -1 , with a similarity matrix U that can be 
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found analytically (but its expression is too lengthy to be 
reported here), and A is the diagonal form of M contain¬ 
ing its eigenvalues Ai = — 2y s , and A 2,3 = — 2y s =F2x|a p |. 
We now solve for the vector v t > ( t ), to find 


3 

Vf{t) = UeW-^U-'M*) = 

n= 1 

(13) 

where we have defined the initial condition vector 

/ Tr s {a|a s A(i')} \ 

UA(t') = v t ’(t') = Tr s {a^(^/)} (i4) 

V Tr s {ai 2 m} ) 


and the matrices M n = [/II n C/ _1 , where II n is a projec¬ 
tor in the n’th “direction”, that is, a matrix with zeros 
everywhere but in element (n, n) which is one. 

Note that for the limit lim^oo v t > (t) to be uniquely 
defined, and therefore for C s to be well-behaved, all the 
eigenvalues of M must satisfy Re{A n } < 0, which in turn 
leads to a requirement for the displacement yoy, < j s . 
This requirement is indeed fulfilled by choosing the mean- 
field displacement as mentioned above, see Sec. |mAl In 
contrast, taking the classical solution as the displacement 
would lead to an ill-behaved C s above and at the classical 
threshold point, that is, for a > 1 . 

Coming back to the correlation functions in eq. 0, 
the general solution (13) allows us to write them all as 


d s (t,t f ) = V]e An(t t) d s , n {t'), 


n= 1 


(15) 


with d s , n {t) = [M n Uj 4 (t )] 3 (the subscript denoting the 
third vector component), where d s denotes any of the cor¬ 
relation functions {d+, d+ ,d~,d~}, for which A is taken, 
respectively, as {5a) 2 p s , p s 5a) 2 , 5a 2 p s , p s 5a 2 }. Let us em¬ 
phasize that, just as with the pump mode, we have been 
able to write all the correlation functions of the signal 
mode into a form where the t dependence only enters in 
simple exponential factors. 

Finally, let us show how this form for the correla¬ 
tion functions allows us to turn the c-MoP equations, 
which are coupled integro-differential equations, into cou¬ 
pled ordinary differential equations. For this aim, let us 
rewrite eqs. © and © as 

Ps(t) = c s p s {t) + | [al 2 (a p )(t) - a 2 s {a p )*(t),Ps(t )] (16) 

+ (f ) 2 {[^,M«)]+h.c.}, 


Pp(t) — £pPp(t) + 2 [ a p( a s) {0 a p( a s)(t) ) Pp(t)\ (17) 


(!) : 


%> 5 ^ ^ ^p,n(^) 


n= 1 


H.c. 


where we have defined the operators 

h s (t) = f dt'e^-^lCsfat^psit'), 

J o 

hp, n (t) = [ dt'e^-^Kp^fat^ppit'), 

Jo 

with the superoperator JC p ^ n defined as IC P in eq. 
but with the correlation functions d s?n (t) instead of d s 
Using their definition, and the solutions found for the 
correlation functions, eqs. (11) and ( p~ 5 ] ) , their evolution 
equations are found to be 


(18) 


0, 

sW- 


d t h s (t ) = (-7 P + C 8 )h 8 {t) +K. 8 {t,t)p 8 (t), (19) 

d t h Pi n(t) = (A n + C p )h P: n(t) + JC Pi n(t, t)p p (t). ( 20 ) 

Together with eqs. ( fl 6 | ) and ( P0 , these form a closed set 
of coupled nonlinear ordinary differential equations for 
the reduced states p s and and the traceless operators 
h s and {h p ^ n } n= 1 ^, 3 . These are the equations that we 
analyze in the remainder of the paper. 

Overall we have shown for the example of the DOPO 
that it is indeed possible to rewrite the integro-differential 
c-MoP equations into a set of ordinary differential equa¬ 
tions. The steps presented here are quite general and 
can be pursued for all c-MoP equations describing any 
physical system. The complexity of the resulting set of 
coupled equations will depend on the complexity of the 
subparts of the full quantum system, here given by the 
complexity of C p and C s . 

Finally, we remark that the c-MoP equations preserve 
the trace and the hermiticity but they do not guarantee 
for the positivity of the density matrix. Such an issue 
is not unusual for projection operator theories, in fact, 
the same conditions can be found in the well established 
Redfield equations [SOI HI]. Obviously whenever the c- 
MoP equations provide a good approximation, they will 
yield a positive density matrix. Hence the positivity of 
the eigenvalues can be used as a consistency test for the 
approximation. 


C. The adiabatic and the diabiatic limit 

The two dissipation rates and set a time scale 
on which the pump and the signal, respectively, relax 
to the steady state of their unperturbed Liouvillians C p 
and C s . In standard open system theory one relies on a 
separation of time scales between the system dynamics 
and the environment correlations. A similar reasoning 
is applied in adiabatic elimination approaches, where for 
the DOPO one relies on a time scale separation between 
signal and pump. The c-MoP theory can, in fact, be 
understood as a generalization of adiabatic elimination 
procedures where one has to consider the back-action of 
the “system” onto the “environment”. We will now show 
that the effective equations for the reduced state of the 
signal known in the adiabatic m and the diabatic |48] 
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limit can, indeed, be obtained as limiting cases of the 
c-MoP equations. 

The adiabatic limit in which the time scale of the pump 
mode is much faster than the time scale of the signal 
mode is defined such that Xp/ls —» oo while XsXp is kept 
finite. The diabatic limit describes the opposite scenario 
where Xp/ls —>> 0. We proceed by comparing the Born 
terms with the free evolution operators C p and £ s , for 
which we consider the scaling of h s /Xs and h p , n / T^which 
can be obtained by simple inspection of eqs. |l9| ) and 
( [ 2 Q| divided by Xs and 7 P , respectively. 

In the adiabatic limit, we infer from eq. p 0 |)/ 7 P that 
h p ,n{t)/ 7 »_= 0 f° r all n and ^ > 0- Introducing this result 
into eq. ( |17| ), we see that the state of the pump will be 
coherent with a field amplitude obeying the equation of 
motion 


dt( a p) — e p 7 p(( a p) T ®p) 2 ^ a s)' (21) 

On the other hand, eq. fiity/xs leads to h s (t)/Xs = 
/C s (t, t)/XsXp = p s (t)5a\ 2 (t)JxsXpi where we have used 
eqs. (p| and ( 11 ) and the fact that when the pump is in a 
coherent state all the expectation values in eq. © can¬ 
cel. Introducing this result into eq. together with 

the steady-state solution of eq. ( 21 ) for (a p ), we end up 
with the effective master equation of the signal mode in 
the adiabatic limit 


7 S 1 dtP s = % [< 


,t2 


at > Ps] + ~^'D a 2p s + V as p s , (22) 


where a = e P x/Xpls is an injection parameter corre¬ 
sponding to a coherent exchange of excitations between 
the signal and pump modes, while g 2 = X 2 /Xpls ac¬ 
counts for signal photon pairs that are lost to the strongly 
damped pump mode. Equation ( [22]) h as been extensively 
studied in the literature mmm- can be derived 
via standard adiabatic elimination which in the language 
of projection operator theory uses a time-independent 
projection superoperator V&d projecting out the coher¬ 
ent laser field ns- Its action on the full density matrix 
is given by V&d pit) = |a)(a| (8) p s (t), where \a) is a co¬ 
herent state with a = e p /x P - The fast exponential decay 


e t ) of the pump correlation functions allows in this 
case for a Markovian approximation in the Born terms, 
that is f* dt’e c °JC s (t,t')p a (t') « K. s {t,t)Ps{t)hp- 
Let us now analyze the c-MoP equations in the dia¬ 
batic limit. In this case, eq. (19^ /x s provides us with 
h s (t)/x s = 0 , which introduced in eq. (16) leads to an 
effective master equation 


Ps(t) = | [( a p )a \ 2 - (a p )*a 2 s ,p s (t)\ +V as p s (t). (23) 


for the signal state. The pump state only enters this 
equation trough the amplitude (a p ) which obeys eq. ( 21 ) 
since h p , n is traceless. Noting that this equation is equiv¬ 


alent to eq. ( 10 j), we conclude that the diabatic limit re¬ 
duces the ful 


1 c-MoP equations to the mean-field equa¬ 


tions. 


We emphasize that within these limits both eqs. (10) 
and ( 22 ) become exact. We have thus shown that the c- 
MoP theory provides us with exact equations of motion 
in the limits Xp/ls oc (adiabatic) and Xp/ls 0 (di¬ 
abatic) where it therefore becomes equivalent with well 
established theories maun. In the remainder of the pa¬ 
per we will step beyond these cases in which time-scale 
separation is present and use the c-MoP theory to access 
the signal state in the 7 p « 7 s scenario. 


IV. ACCURACY TESTS AND FULL QUANTUM 
STATES OF THE SIGNAL MODE 


In the previous section we have shown how to deal 
with the non-Markovian structure of the c-MoP equa¬ 
tions. The only remaining difficulty is given by the quar- 
tic structure of the effective equations of motion derived 
for the signal mode, eqs. (16) and (19). In this section we 
will treat the problem numerically in the Fock state basis 
by introducing a truncation D s for the Hilbertspace H s 
of the signal, where D s is chosen such that the results 
for the observables we are interested in converge up to 
some desired accuracy. Thus, the reduced state p s and 
the operator h s (t) will be D s x D s dimensional matri¬ 
ces. Instead of treating the pump mode in an analogous 
manner, we exploit the fact that the c-MoP equations 
of the pump mode © and ( |20| ) are quadratic in the 
bosonic operators. As a consequence we are able to de¬ 
scribe the pump state by a set of closed equations for five 
variables only, the mode amplitude (a p ) plus the fluctu¬ 
ations ( a p Sa p ) and (a^ p Sa p ) (note that the first two are 
complex variables). At the end, we are thus effectively 
left with two coupled differential equations for the ma¬ 
trices p s (t) and h s (t ), with the pump equations solved 
either in parallel numerically or analytically as a func¬ 
tion of signal observables. 

In this section, we will compare the steady states of 
the classical theory from eq. the steady states of the 
mean-field equations ( JI(j|), and the steady states of the c- 
MoP equations 0 and ([ 5 ]) . In order to show the accuracy 
of the c-MoP equations we will also determine the steady 
state of the full Liouville-von Neumann equation ([l]) in 
parameter regimes where it is numerically tractable. This 
numerical simulation is done as follows: first, we elimi¬ 
nate the large coherent background of the laser drive from 
the Liouvillian C by writing a p = a p + Sa p , where a p i s 


taken to be the classical steady-state solution of eqs. 
then, we use the superspace formalism, where the steady- 
state operator p ss and the Liouville superoperator C are 
represented, respectively, by a vector p ss and a matrix L, 
and p ss can be found as the eigenvector with zero eigen¬ 
value of L [52, 53] . As the dimension of the matrix L is 
(D P x D s ) 2 , with D p denoting the pump mode’s Hilbert 
space dimension, this exact simulation is limited to small 
photon numbers. 

In all the simulations we consider cases without time- 
scale separation between the two modes and rescale all 















FIG. 2. Accuracy tests of the c-MoP theory for steady-state expectation values as a function of the injection parameter a. 
In all plots we set j p = j 8 = 1. The cases x — 1 and X — 0.1 are considered in (a)-(c) and (<f)-(/), respectively. The rescaled 
pump amplitude x( a p) 1S shown in ( a ) and (d); ( b ) and (e) show the signal photon number (aja s ); finally, (c) and (/) show 
the g^ function of the signal, which is equal to 1 for a coherent state (or a balanced mixture of them). The gray thin solid 
lines show the classical prediction from eqs. showing that the classical threshold where the signal field is switc hed on lies 
at a = 1. The blue solid curves represent the results obtained from the numerical solution of the c-MoP equations (16), |l7| ), 
(19), and (20). The red stars show the result obtained from the full master equation 0 up to injection parameters a where 
the numerics are tractable for us. Finally, the black dashed curves represent the mean-field theory, see eq. (10). Apart from 
the classical solution, all theories conserve the Z 2 symmetry, i.e. (a s ) — 0. 


units to the dissipation rates, i.e. we put 7 P = = 1. 

The only remaining parameters are the nonlinear cou¬ 
pling % and the injection parameter a — e P x- 

In Fig.[2]we present results in parameter regimes where 
the full DOPO equation can be solved numerically. 
In Figs. [2Ja) — (c) and [2Jd) — (/) we show different 
steady-state observables for x = 1 and x = 0-1? respec¬ 
tively. It can be appreciated how the c-MoP results (blue 
solid line) coincide almost perfectly with the numerical 
results from the full master equation (red stars). The 
observables that we show are the pump mode’s ampli¬ 
tude (a p ) in Figs. [2fa) andl^d), the signal photon num¬ 
ber (a|a s ) in Figs7l2[6) and]2[e), and the g^ function 
^^(O) = (a| 2 a 2 )/(a|a s ) 2 of the signal in Figs, [^c) and 
[2J/). We also compare with the mean-field predictions 
of eqs. ( flQ| ) (black dashed line), which in this context 
should be understood as the c-MoP theory up to first 
order, and with the classical steady-state solutions (gray 
thin solid line) given after eq. ©• Let us remark that 
despite the nonlinear nature of the mean-field and the 
c-MoP equations, we only find one physical solution for 
each of them. 

All four theories agree quite well far below the crit¬ 
ical point a = 1 as the states of the signal and pump 
modes are close to vacuum and a coherent state induced 
by the external laser drive, respectively. Far above the 
threshold point, where the classical theory is expected 
to be approximately valid, we find that both the c-MoP 


predictions and the exact numerics agree well with the 
classical solutions for all observables, but with the funda¬ 
mental difference that the classical theory breaks the Z 2 
symmetry, while c-MoP and the exact solution preserve 
it. The mean-field solution, on the other hand, fails to 
describe the state of the signal above threshold as can 
be appreciated from the g^ function in Figs, [^c) and 
0 /)- As expected, mean-field theory and the classical 
theory break down in the vicinity of the threshold point. 
Remarkably, this is not true for c-MoP which appears to 
give quasi exact results for all values of <7, even in cases 
where the interaction rate x is comparable to all other 
system parameters. 

For the experimentally relevant scenario with x 1, 
the Hilbert space dimension needs to be so large that 
we are not able to find the numerical solution of the 
full master equation 0 for injection parameters close 
to (or above) threshold. However, we can compare the 
c-MoP predictions (red stars), see Fig. with the per¬ 
turbative approach which Drummond et al. (dark yellow 
dot-dashed line) developed in the vicinities of the critical 
point, by making a consistent multiple-scale expansion 
of the system’s stochastic variables within the positive 
P representation [3H [35] . This procedure has the virtue 
of being valid for any values of and y s , and close to 
threshold, concretely for \a — 1| < x/ s-, it is ex¬ 
pected to be quasi-exact. As shown in Figs. |4^a) and 
m we find perfect agreement between this approach 
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FIG. 3. Wigner functions of the c-MoP density matrix for the signal mode without (a) and with (b) the Gaussian state 
approximation for = 1, % = 0.1 and for different values of a. In the absence of injection, a = 0, the signal state is in 

vacuum. Upon approaching the threshold, it becomes squeezed, with the highest squeezing levels obtained around a = 1. Above 
threshold two symmetric peaks appear and the squeezing reaches some asymptotic value as we move away from threshold. Note 
how above threshold the state can be approximated by a balanced mixture of two symmetry breaking states. Indeed, let us 
remark that while for a < 1 we are plotting the unique solution that appears when applying the Gaussian state approximation 
onto the c-MoP equations (which we have called below threshold solution in the text), for a > 1 we have chosen to plot the 
Wigner function corresponding to a balanced mixture of the two above threshold symmetry breaking solutions with opposite 
phase which coexist with the non-symmetry breaking Gaussian solution. 


and the c-MoP theory for x = 0.01. 

Overall, we have indeed shown the drastic impact of 
the Born terms, which do not only lead to a quantitative 
improvement as compared to the classical theory or to 
mean-field, but to a qualitatively different state of the 
signal mode. The classical theory predicts a coherent 
state, while the mean-field theory, i.e. the c-MoP theory 
up to first order, predicts a Gaussian state of the sig¬ 
nal centered around (a s ) = 0 [43]. The c-MoP theory 
including the born terms, hence including signal-pump 
correlations within a projection operator based theory, 
is capable of finding the full quantum state of the signal 
which is neither coherent nor Gaussian as shown through 
the function in Figs. life) and Iff). 

In order to illustrate the full quantum state, we plot 
the Wigner function W(x s ,p s ) of the signal density ma¬ 
trix obtained from the c-MoP equations in Fig. [3Ja) for 
X = 0.1 and different values of a. Let us remark that 
in our case in which the Wigner function is positive ev¬ 
erywhere in the phase-space formed by the quadratures 
x s = a| + a s and p s = i(a\ — a s ), it can be simply 
interpreted as the joint probability distribution describ¬ 
ing the statistics of measurements of these observables 
m nai ng. From a computational point of view, we 
evaluate it from the steady-state density matrix follow¬ 
ing the method detailed in [54] . Far below threshold, the 
Wigner function shows a perfect vacuum for the signal 


state, see top panel of Fig. [3Ja) for a = 0 as a refer¬ 
ence. As we cross through the critical point, two sig¬ 
nificant effects take place. First, approaching the thresh¬ 
old we find the well-known quadrature-noise reduction or 
squeezing HJ3IE3G3I, which is highest around the crit¬ 
ical point a — 1 [35], and reaches its asymptotic value 
(Spl) = ( 7 s + 7 p)/( 27 s +7 P ) for a —)> oo, with correspond¬ 
ing antisqueezing (Sx 2 ) = 1 +7 S /7 P [43]. Second, as we 
cross the threshold we appreciate how the state devel¬ 
ops two peaks centered (asymptotically) at the quadra¬ 
ture values predicted by the classical solution. Hence, 
even though the true quantum state never breaks the Z 2 
symmetry, it does so in two qualitatively different ways 
depending on whether we are below or above threshold. 

In retrospect, we see that the symmetry-breaking 
states predicted above threshold by the standard lin¬ 
earization approach correspond each to one of the two 
distinct peaks appearing in the exact state. Far above 
threshold a 1 the two peaks have zero overlap and 
such states provide reasonable predictions for all observ¬ 
ables which are not sensitive to symmetry breaking, that 
is, all observables containing even numbers of signal field 
operators. Of course, such a deficit can be corrected 
by simply using a balanced mixture of the symmetry¬ 
breaking states [43]; this construction will guide us in the 
next section, where we will perform a Gaussian approxi¬ 
mation which necessarily breaks the Z 2 symmetry. It is 
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then close to the critical point where both linearization 
and mean-field approaches fail, whereas c-MoP provides 
an accurate description of the quantum state. 

Let us remark that we have compared the Wigner 
function obtained from the c-MoP theory with the re¬ 
duced signal states obtained from the full master equa¬ 
tion, which was only possible for a < 1.2, and found 
very good agreement, the differences being completely 
unnoticeable to the naked eye. We emphasize again that, 
with the numerical solution of the c-MoP equations we 
are able to find the full reduced density matrices of the 
modes away from the adiabatic limit. This is in con¬ 
trast to other approaches such as stochastic simulations 
j34b37] or the Keldysh formalism [3814401 which are nat¬ 
urally design to provide expectation values of the system 
operators. 


V. GAUSSIAN STATE APPROXIMATION 
WITHIN THE C-MOP THEORY 


Despite the fact that the complexity of solving the c- 
MoP equations fully numerically scales in a more favor¬ 
able way than the numerical complexity of the full master 
equation, it still requires to integrate a number of differ¬ 
ential equations that scales quadratically with the dimen¬ 
sion of the truncated Hilbert space for the signal field. 
Therefore, it is very desirable to find an effective descrip¬ 
tion of the underlying theory which is numerically more 
efficient and can thus cover the whole parameter space. 
In the remainder of this section, we implement such an 
idea by applying a Gaussian state approximation (GSA) 
consistent with the c-MoP equations 0 and 0. 

Another great advantage of a Gaussian theory, apart 
from reaching the whole parameter space, is the efficiency 
in the evaluation of both steady states and dynamical 
quantities such as two-time correlation functions. The 
disadvantage of a Gaussian theory, however, is the lack 
of quantitative accuracy especially in the vicinity of the 
critical point. Nonetheless, as we show in the follow¬ 
ing, a Gaussian theory consistent with the c-MoP equa¬ 
tions offers better quantitative accuracy than any of the 
previously developed Gaussian methods, particularly lin¬ 
earization around the classical solution or the recently- 
developed self-consistent linearization [43] . 

The general procedure for finding a GSA for the state 
of a certain bosonic master equation is very simple. In a 
first step, we write the bosonic operators as aj = ctj+Saj, 
with Oij — {a,j), such that (Saj) = 0. In the next step 
we find the evolution equation for the first and second 
moments, which will depend on higher-order moments in 
general. Thus, in the final step we assume the state to 
be Gaussian at all times, so that all higher order mo¬ 
ments factorize into products of first and second order 
moments [H]|43]; in particular, we will encounter third 
order moments such as e.g. ( Sal 2 Sa s ) which vanish iden¬ 
tically within the GSA, and forth order moments which 


factorize according to, e.g. 

(5«t 4 ) * 3(<5aJ 2 ) 2 , 

(Sal 3 Sa s ) *3{6at 2 ){6at6a s ), (24) 

(Sal 2 5a 2 s ) « (Sal 2 ) (5^) + 2(da\da s ) 2 . 


After this final step, we are then left with a closed set of 
nonlinear equations for the amplitudes ay and the second 
order moments of the fluctuations Saj that have to be 
solved self-consistently. 

The standard linearization theory can be understood 
as a GSA on the full master equation, but with the ex¬ 
ception that the amplitudes ay are not determined self- 
consistently, but are obtained from the classical theory. 
As shown by the gray thin solid line in Fig. [4 Jb) the 
complete suppression of quantum fluctuations when de¬ 
termining these amplitudes leads to unphysical results at 
the threshold point in the DOPO. 

The self-consistent linearization method , as it is coined 
in Reference [43], goes one step beyond standard lin¬ 
earization by consistently finding the amplitudes ay from 
the GSA still applied to the full master equation. Due 
to the nonlinear nature of the resulting equations of mo¬ 
tion one can find several solutions in a given point of 
parameter space. However, it was shown that at the end 
only two types of solutions were physical, qualitatively 
similar to the solutions found from standard lineariza¬ 
tion, but quantitatively regularized in such a way that 
the unphysical results of the latter disappear. In partic¬ 
ular, a below threshold (BT) solution was found, which 
does not break the Z 2 symmetry, i.e. a s = 0, but in 
contrast to the classical theory exists for all values of the 
injection parameter, not only for a < 1. We also found 
two above threshold (AT) solutions with opposite phase 
which break the symmetry, that is, (a s ) = ±\ot s \ 7 ^ 0, but 
appear only above a certain injection parameter a > 1 
which is larger than the classical threshold value. Inter¬ 
estingly, we point out that the BT solution found through 
this self-consistent linearization is exactly equivalent to 
the mean-field theory introduced in Section [III A| 

Motivated by these findings, we apply a GSA to the 
c-MoP equations. Concretely, we calculate all first and 
second order moments of the pump and signal fluctua¬ 
tions from the c-MoP eqs. (16]), ©, and < [2Q| ), and 

apply the factorization of higher order moments as ex¬ 
plained above. In strong contrast to the GSA on the full 
master equation, we do not need to assume a Gaussian 
state for the full state p but only for the reduced state of 
the signal p s . Hence, we expect similar qualitative results 
but with a higher quantitative accuracy. 

Indeed, this is what we find and illustrate in Fig. [4] for 
Ip = 7s = 1 and % = 0.01. We plot steady state expec¬ 
tation values for the pump amplitude x( a p) 1 n Fig. 0a), 
the signal photon number (a|a s ) = (5a|<5a s ) + | a Jr in 
Fig. ]4|6), and the g^ function of the signal in Fig. ]4[c), 
all as a function of the injection parameter a. The blue 
solid line in Fig. [4] shows the below threshold solution of 
the GSA on the c-MoP equation, while the green dashed 
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FIG. 4. Accuracy tests of the c-MoP theory with and without the Gaussian state approximation for steady-state observables 
as a function of the injection a. In all plots we set 7 P = = 1 and x — 0-01- As in fig. [2] we show the rescaled pump 

field amplitude x( a p)> the signal photon number (oJo s ), and the g ^ function of the signal, in (a), (6), and (c), respectively. 
The red stars show the results obtained from the c-MoP eqs. ( [16] ), ( |19| , and ( |20| ), up to injection parameters a where 

the numerics are tractable. For comparison, the quasi exact method of Drummond and collaborators El EH is shown as a 
dark yellow dot-dashed line. The blue solid and the green dashed curves represent the below and above threshold solutions, 
respectively, obtained from a Gauss ian state approximation on the c-MoP equations. The black thin dotted curve displays the 
results of mean-field theory, see eq. ( [To] ) , which in this context can be understood as the below threshold solution of a Gaussian 
state approximation on the full master equation 0- Finally, the gray thin solid lines represent the prediction of the standard 
linearization theory in (a) and (6), and the coherent-state prediction g^ = 1 of the classical equations ([ 2 ]) in (c). 


line illustrates the above threshold solution. The latter 
fulfills (Sa\Sa s ) <C |<a s | 2 and is therefore more likely to 
provide physically consistent results than the BT solution 
whenever they coexist. In Fig.Qc) we show how the AT 
solution indeed gives the correct value for the g^ func¬ 
tion, what indicates that each of the AT solutions corre¬ 
sponds to one of the peaks of the Wigner function, see 
Fig. J3^a), and considers Gaussian fluctuations around it. 
In order to illustrate this point even further, we show in 
Fig. [3Jb) the Wigner function [12, T3j 09] corresponding 
to the GSA on the c-MoP equations (as explained in the 
previous section, above threshold we take the balanced 
mixture of the two symmetry breaking solutions, such 
that the resulting state preserves the Z 2 symmetry). 

Importantly, there is an increased quantitative accu¬ 
racy of the BT solution obtained from the c-MoP the¬ 
ory as compared with the mean-field theory (or the self- 
consistent linearization), see Figs. [4Ja) and [4^6), for pa¬ 
rameters below and especially at the classical threshold 
point. As mentioned in Sec. ( |lV| ) we test the accuracy of 
our method by comparing with the quasi exact method of 
Drummond and collaborators pi [Mi illustrated by the 
dark yellow dot-dashed line in Figs. [4ja) and [4^6). This 
increase in accuracy can be attributed to the born terms, 
since the mean-field equations can be understood on the 
one hand as the first order approximation of the c-MoP 
theory, and on the other hand as the below threshold 
solution of the GSA on the full master equation of the 
DOPO. 

To summarize this section, we have shown that the c- 
MoP equations also provide a highly accurate Gaussian 
theory which is still as effective as every other linearized 
theory but, in contrast, it takes significant signal-pump 
correlations into account. This is relevant because, as 
stated above, a Gaussian theory has the virtue that both 
steady-state as well as dynamical quantities such as two- 
time correlation functions can be found efficiently for any 


time and set of parameters. To emphasize this practical 
aspect of the GSA, we will show in Sec. |VI| that the level 
of accuracy that we have found here in the evaluation 
of the steady states is also present in the transient time 
evolution. 


VI. DYNAMICS 

So far we have only presented steady state quantities 
for the various methods of our interest. In this section 
we will briefly elaborate on the possibility to simulate 
dynamical evolution as well. The steady state of the full 
master equation 0 can be understood as an eigenvector 
corresponding to the zero eigenvalue of the Liouvillian C 
such that Cp ss = 0. The formal solution for the time 
evolving state which can be written as p(t) = e Ct p( 0), on 
the other hand, involves all eigenvalues of the Liouvillian. 
Hence, it is a priori not clear whether a given approxi¬ 
mate method used for the evaluation of the steady state 
of C will provide the same degree of accuracy when used 
for transient time evolution. 

In order to investigate this open issue we simulate the 
time dynamics of the various approximate methods that 
we have introduced, and compare their results with an 
exact simulation of the full master equation ([l]) in regions 
of the parameter space where it is numerically tractable. 
We analyze a situation in which the input laser drives 
the system from the initial vacuum to its steady state. 
Figure [5] shows the signal photon number as a func¬ 
tion of time at the classical threshold point a = 1, for 
Ip = 7s = 1, an d for x = 0.1 in Figure j5^a) and x = 0-05 
in Figure [5^6). The red stars in Fig. |5[a) illustrate the 
result obtained from the numerical simulation of the full 
master equation 0, while the red line in Fig. §6) il- 
lustrates the steady state value of the observables only, 
since the small value of x prevented us from being able to 
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FIG. 5. Accuracy tests of the c-MoP theory for transient time evolution. The initial state is chosen to be the vacuum. We set 
7 p = 7s = 1, investigate the classical threshold point a = 1, and choose x — 0-1 in ( a ) and X — 0-05 in ( b ). The plots display 
the signal photon number as a function of time in units of the dissipation rates, while the insets show the g ^ function of the 
signal mode. The red stars in (a) show the result obtained from the numerical simulation of the full master equation |l]), while 
the red line in ( b ) indicates its steady-state values only. The blue solid curves represent the results obtained from the numerical 
integration of the c-MoP eqs. ( [16] ), <Ezf |l9] ), and ( [20| . Finally, the green dashed and black dotted lines represent the time 
evolution obtained from a Gaussian state approximation on the c-MoP equations and the full master equation (mean-field), 
respectively. 


simulate the whole dynamics in this case. On the other 
hand, the blue solid curves represent the results obtained 
from the numerical integration of the c-MoP equations 
as explained in Section [IV] Finally, the green dashed and 
black dotted lines represent the time evolution obtained 
from a GSA on the c-MoP equations and the full master 
equation, respectively. 

Remarkably, Fig.|5|a) shows that the level of accuracy 
found dynamically for the various approximations is sim¬ 
ilar to the ones that we already encountered when evalu¬ 
ating steady-state quantities. In particular, it is apparent 
that, at any point in time, the GSA on the full master 
equation is less accurate than the GSA on the c-MoP 
equations, which in turn does not have the remarkable 
level of accuracy shown by the full c-MoP numerical sim¬ 
ulation, which almost coincides with the numerics of the 
full master equation at all times. It is important to note 
that the evolution of the g^ function shown in the inset 
of Fig. [5|a) suggests that, indeed, the c-MoP equations 
are able to map the full quantum state of the signal in 
the course of time. 

A numerical simulation for the parameter set chosen in 
Fig. [5Jb) demands minimal Hilbert space dimensions of 
dim Up = 6 and dimH s = 120 in order to reach conver¬ 
gence up to an accuracy of 10“ 2 for the relevant observ¬ 
ables. Thus, while the c-MoP approach requires a sim¬ 
ulation of a set of 28 811 coupled nonlinear differential 
equations, in the case of the full master equation one has 
to integrate 518 400 coupled linear differential equations, 
which has precluded us from being able to simulate the 
dynamics from it. Therefore we only show steady-state 
observables of the full master equation for these case. 

Figs.[5ja) and[5|6) further illustrate the scaling of var¬ 
ious quantities with the nonlinear coupling x at the criti¬ 
cal point. In particular, note how both the signal photon 


number and the time that the system needs to reach the 
steady state double when x is reduced by half. The latter 
is known in the literature as critical slowing down m, 
and just as the signal photon number, it was predicted to 
scale with x" 1 HHIMIEH], in agreement with our c-MoP 
simulation. Hence, we can appreciate the practical use 
of a Gaussian theory by considering that to simulate an 
experimentally relevant scenario where x "C 1? dynam¬ 
ical quantities would require extremely long simulation 
times, which, as explained before, can be efficiently han¬ 
dled with a GSA on the c-MoP theory, but not by its full 
numerical simulation. As an example, we have checked 
that for x = 0.01 a GSA on the c-MoP equations re¬ 
quires a normalized time of approximately 300 to reach 
the steady state, again in agreement with the x -1 scal¬ 
ing, as it can be appreciated in Fig. [5ja) that such time 
is about 10 times smaller for x = 0.1. 

VII. CONCLUSIONS AND OUTLOOK 

In conclusion we have exemplified the applicability of 
the self-consistent projection operator theory to nonlin¬ 
ear quantum optical systems on the case study of the de¬ 
generate optical parametric oscillator. Our theory gener¬ 
alizes mean-field approaches and in particular adiabatic 
elimination methods to settings without time-scale sepa¬ 
ration. The effective master equations can be solved effi¬ 
ciently despite their non-Markovian structure. We have 
demonstrated the high degree of accuracy of our method 
and revealed its capability to determine the exact quan¬ 
tum states below, at, and above the classical threshold 
for both stationary states and dynamical evolution. 

In addition, we developed a linearized theory consis¬ 
tent with the self-consistent Mori projector equations and 
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showed its accuracy far beyond any known linearized the¬ 
ories. We expect our Gaussian method to be particularly 
useful in the context of hybrid systems such as optome¬ 
chanical parametric oscillators [25, 26], where fields of 
quantum nature with no coherent background are cou¬ 
pled to mechanical elements. Some intriguing tasks for 
future research would include applying the c-MoP ap¬ 
proach to investigate dynamical questions, e.g. investi¬ 
gate tunneling times between the two symmetry breaking 
states in parameter regimes away from the adiabatic limit 
[52], simulate quantum quenches in a driven-dissipative 
scenario m, and investigate the effect of small symmetry 
breaking perturbations on both the dynamics and steady 


states. 
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